// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012-2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
// Copyright (C) 2012-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_SPARSE_QR_H
#define EIGEN_SPARSE_QR_H

namespace Eigen {

template<typename MatrixType, typename OrderingType>
class SparseQR;
template<typename SparseQRType>
struct SparseQRMatrixQReturnType;
template<typename SparseQRType>
struct SparseQRMatrixQTransposeReturnType;
template<typename SparseQRType, typename Derived>
struct SparseQR_QProduct;
namespace internal {
template<typename SparseQRType>
struct traits<SparseQRMatrixQReturnType<SparseQRType>>
{
	typedef typename SparseQRType::MatrixType ReturnType;
	typedef typename ReturnType::StorageIndex StorageIndex;
	typedef typename ReturnType::StorageKind StorageKind;
	enum
	{
		RowsAtCompileTime = Dynamic,
		ColsAtCompileTime = Dynamic
	};
};
template<typename SparseQRType>
struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType>>
{
	typedef typename SparseQRType::MatrixType ReturnType;
};
template<typename SparseQRType, typename Derived>
struct traits<SparseQR_QProduct<SparseQRType, Derived>>
{
	typedef typename Derived::PlainObject ReturnType;
};
} // End namespace internal

/**
 * \ingroup SparseQR_Module
 * \class SparseQR
 * \brief Sparse left-looking QR factorization with numerical column pivoting
 *
 * This class implements a left-looking QR decomposition of sparse matrices
 * with numerical column pivoting.
 * When a column has a norm less than a given tolerance
 * it is implicitly permuted to the end. The QR factorization thus obtained is
 * given by A*P = Q*R where R is upper triangular or trapezoidal.
 *
 * P is the column permutation which is the product of the fill-reducing and the
 * numerical permutations. Use colsPermutation() to get it.
 *
 * Q is the orthogonal matrix represented as products of Householder reflectors.
 * Use matrixQ() to get an expression and matrixQ().adjoint() to get the adjoint.
 * You can then apply it to a vector.
 *
 * R is the sparse triangular or trapezoidal matrix. The later occurs when A is rank-deficient.
 * matrixR().topLeftCorner(rank(), rank()) always returns a triangular factor of full rank.
 *
 * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<>
 * \tparam _OrderingType The fill-reducing ordering method. See the \link OrderingMethods_Module
 *  OrderingMethods \endlink module for the list of built-in and external ordering methods.
 *
 * \implsparsesolverconcept
 *
 * The numerical pivoting strategy and default threshold are the same as in SuiteSparse QR, and
 * detailed in the following paper:
 * <i>
 * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
 * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011.
 * </i>
 * Even though it is qualified as "rank-revealing", this strategy might fail for some
 * rank deficient problems. When this class is used to solve linear or least-square problems
 * it is thus strongly recommended to check the accuracy of the computed solution. If it
 * failed, it usually helps to increase the threshold with setPivotThreshold.
 *
 * \warning The input sparse matrix A must be in compressed mode (see SparseMatrix::makeCompressed()).
 * \warning For complex matrices matrixQ().transpose() will actually return the adjoint matrix.
 *
 */
template<typename _MatrixType, typename _OrderingType>
class SparseQR : public SparseSolverBase<SparseQR<_MatrixType, _OrderingType>>
{
  protected:
	typedef SparseSolverBase<SparseQR<_MatrixType, _OrderingType>> Base;
	using Base::m_isInitialized;

  public:
	using Base::_solve_impl;
	typedef _MatrixType MatrixType;
	typedef _OrderingType OrderingType;
	typedef typename MatrixType::Scalar Scalar;
	typedef typename MatrixType::RealScalar RealScalar;
	typedef typename MatrixType::StorageIndex StorageIndex;
	typedef SparseMatrix<Scalar, ColMajor, StorageIndex> QRMatrixType;
	typedef Matrix<StorageIndex, Dynamic, 1> IndexVector;
	typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
	typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;

	enum
	{
		ColsAtCompileTime = MatrixType::ColsAtCompileTime,
		MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
	};

  public:
	SparseQR()
		: m_analysisIsok(false)
		, m_lastError("")
		, m_useDefaultThreshold(true)
		, m_isQSorted(false)
		, m_isEtreeOk(false)
	{
	}

	/** Construct a QR factorization of the matrix \a mat.
	 *
	 * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()).
	 *
	 * \sa compute()
	 */
	explicit SparseQR(const MatrixType& mat)
		: m_analysisIsok(false)
		, m_lastError("")
		, m_useDefaultThreshold(true)
		, m_isQSorted(false)
		, m_isEtreeOk(false)
	{
		compute(mat);
	}

	/** Computes the QR factorization of the sparse matrix \a mat.
	 *
	 * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()).
	 *
	 * \sa analyzePattern(), factorize()
	 */
	void compute(const MatrixType& mat)
	{
		analyzePattern(mat);
		factorize(mat);
	}
	void analyzePattern(const MatrixType& mat);
	void factorize(const MatrixType& mat);

	/** \returns the number of rows of the represented matrix.
	 */
	inline Index rows() const { return m_pmat.rows(); }

	/** \returns the number of columns of the represented matrix.
	 */
	inline Index cols() const { return m_pmat.cols(); }

	/** \returns a const reference to the \b sparse upper triangular matrix R of the QR factorization.
	 * \warning The entries of the returned matrix are not sorted. This means that using it in algorithms
	 *          expecting sorted entries will fail. This include random coefficient accesses (SpaseMatrix::coeff()),
	 *          and coefficient-wise operations. Matrix products and triangular solves are fine though.
	 *
	 * To sort the entries, you can assign it to a row-major matrix, and if a column-major matrix
	 * is required, you can copy it again:
	 * \code
	 * SparseMatrix<double>          R  = qr.matrixR();  // column-major, not sorted!
	 * SparseMatrix<double,RowMajor> Rr = qr.matrixR();  // row-major, sorted
	 * SparseMatrix<double>          Rc = Rr;            // column-major, sorted
	 * \endcode
	 */
	const QRMatrixType& matrixR() const { return m_R; }

	/** \returns the number of non linearly dependent columns as determined by the pivoting threshold.
	 *
	 * \sa setPivotThreshold()
	 */
	Index rank() const
	{
		eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
		return m_nonzeropivots;
	}

	/** \returns an expression of the matrix Q as products of sparse Householder reflectors.
	 * The common usage of this function is to apply it to a dense matrix or vector
	 * \code
	 * VectorXd B1, B2;
	 * // Initialize B1
	 * B2 = matrixQ() * B1;
	 * \endcode
	 *
	 * To get a plain SparseMatrix representation of Q:
	 * \code
	 * SparseMatrix<double> Q;
	 * Q = SparseQR<SparseMatrix<double> >(A).matrixQ();
	 * \endcode
	 * Internally, this call simply performs a sparse product between the matrix Q
	 * and a sparse identity matrix. However, due to the fact that the sparse
	 * reflectors are stored unsorted, two transpositions are needed to sort
	 * them before performing the product.
	 */
	SparseQRMatrixQReturnType<SparseQR> matrixQ() const { return SparseQRMatrixQReturnType<SparseQR>(*this); }

	/** \returns a const reference to the column permutation P that was applied to A such that A*P = Q*R
	 * It is the combination of the fill-in reducing permutation and numerical column pivoting.
	 */
	const PermutationType& colsPermutation() const
	{
		eigen_assert(m_isInitialized && "Decomposition is not initialized.");
		return m_outputPerm_c;
	}

	/** \returns A string describing the type of error.
	 * This method is provided to ease debugging, not to handle errors.
	 */
	std::string lastErrorMessage() const { return m_lastError; }

	/** \internal */
	template<typename Rhs, typename Dest>
	bool _solve_impl(const MatrixBase<Rhs>& B, MatrixBase<Dest>& dest) const
	{
		eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
		eigen_assert(this->rows() == B.rows() &&
					 "SparseQR::solve() : invalid number of rows in the right hand side matrix");

		Index rank = this->rank();

		// Compute Q^* * b;
		typename Dest::PlainObject y, b;
		y = this->matrixQ().adjoint() * B;
		b = y;

		// Solve with the triangular matrix R
		y.resize((std::max<Index>)(cols(), y.rows()), y.cols());
		y.topRows(rank) =
			this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank));
		y.bottomRows(y.rows() - rank).setZero();

		// Apply the column permutation
		if (m_perm_c.size())
			dest = colsPermutation() * y.topRows(cols());
		else
			dest = y.topRows(cols());

		m_info = Success;
		return true;
	}

	/** Sets the threshold that is used to determine linearly dependent columns during the factorization.
	 *
	 * In practice, if during the factorization the norm of the column that has to be eliminated is below
	 * this threshold, then the entire column is treated as zero, and it is moved at the end.
	 */
	void setPivotThreshold(const RealScalar& threshold)
	{
		m_useDefaultThreshold = false;
		m_threshold = threshold;
	}

	/** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
	 *
	 * \sa compute()
	 */
	template<typename Rhs>
	inline const Solve<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const
	{
		eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
		eigen_assert(this->rows() == B.rows() &&
					 "SparseQR::solve() : invalid number of rows in the right hand side matrix");
		return Solve<SparseQR, Rhs>(*this, B.derived());
	}
	template<typename Rhs>
	inline const Solve<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const
	{
		eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
		eigen_assert(this->rows() == B.rows() &&
					 "SparseQR::solve() : invalid number of rows in the right hand side matrix");
		return Solve<SparseQR, Rhs>(*this, B.derived());
	}

	/** \brief Reports whether previous computation was successful.
	 *
	 * \returns \c Success if computation was successful,
	 *          \c NumericalIssue if the QR factorization reports a numerical problem
	 *          \c InvalidInput if the input matrix is invalid
	 *
	 * \sa iparm()
	 */
	ComputationInfo info() const
	{
		eigen_assert(m_isInitialized && "Decomposition is not initialized.");
		return m_info;
	}

	/** \internal */
	inline void _sort_matrix_Q()
	{
		if (this->m_isQSorted)
			return;
		// The matrix Q is sorted during the transposition
		SparseMatrix<Scalar, RowMajor, Index> mQrm(this->m_Q);
		this->m_Q = mQrm;
		this->m_isQSorted = true;
	}

  protected:
	bool m_analysisIsok;
	bool m_factorizationIsok;
	mutable ComputationInfo m_info;
	std::string m_lastError;
	QRMatrixType m_pmat;			// Temporary matrix
	QRMatrixType m_R;				// The triangular factor matrix
	QRMatrixType m_Q;				// The orthogonal reflectors
	ScalarVector m_hcoeffs;			// The Householder coefficients
	PermutationType m_perm_c;		// Fill-reducing  Column  permutation
	PermutationType m_pivotperm;	// The permutation for rank revealing
	PermutationType m_outputPerm_c; // The final column permutation
	RealScalar m_threshold;			// Threshold to determine null Householder reflections
	bool m_useDefaultThreshold;		// Use default threshold
	Index m_nonzeropivots;			// Number of non zero pivots found
	IndexVector m_etree;			// Column elimination tree
	IndexVector m_firstRowElt;		// First element in each row
	bool m_isQSorted;				// whether Q is sorted or not
	bool m_isEtreeOk;				// whether the elimination tree match the initial input matrix

	template<typename, typename>
	friend struct SparseQR_QProduct;
};

/** \brief Preprocessing step of a QR factorization
 *
 * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()).
 *
 * In this step, the fill-reducing permutation is computed and applied to the columns of A
 * and the column elimination tree is computed as well. Only the sparsity pattern of \a mat is exploited.
 *
 * \note In this step it is assumed that there is no empty row in the matrix \a mat.
 */
template<typename MatrixType, typename OrderingType>
void
SparseQR<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
{
	eigen_assert(
		mat.isCompressed() &&
		"SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
	// Copy to a column major matrix if the input is rowmajor
	typename internal::conditional<MatrixType::IsRowMajor, QRMatrixType, const MatrixType&>::type matCpy(mat);
	// Compute the column fill reducing ordering
	OrderingType ord;
	ord(matCpy, m_perm_c);
	Index n = mat.cols();
	Index m = mat.rows();
	Index diagSize = (std::min)(m, n);

	if (!m_perm_c.size()) {
		m_perm_c.resize(n);
		m_perm_c.indices().setLinSpaced(n, 0, StorageIndex(n - 1));
	}

	// Compute the column elimination tree of the permuted matrix
	m_outputPerm_c = m_perm_c.inverse();
	internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
	m_isEtreeOk = true;

	m_R.resize(m, n);
	m_Q.resize(m, diagSize);

	// Allocate space for nonzero elements: rough estimation
	m_R.reserve(2 *
				mat.nonZeros()); // FIXME Get a more accurate estimation through symbolic factorization with the etree
	m_Q.reserve(2 * mat.nonZeros());
	m_hcoeffs.resize(diagSize);
	m_analysisIsok = true;
}

/** \brief Performs the numerical QR factorization of the input matrix
 *
 * The function SparseQR::analyzePattern(const MatrixType&) must have been called beforehand with
 * a matrix having the same sparsity pattern than \a mat.
 *
 * \param mat The sparse column-major matrix
 */
template<typename MatrixType, typename OrderingType>
void
SparseQR<MatrixType, OrderingType>::factorize(const MatrixType& mat)
{
	using std::abs;

	eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
	StorageIndex m = StorageIndex(mat.rows());
	StorageIndex n = StorageIndex(mat.cols());
	StorageIndex diagSize = (std::min)(m, n);
	IndexVector mark((std::max)(m, n));
	mark.setConstant(-1);		  // Record the visited nodes
	IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q
	Index nzcolR, nzcolQ;		  // Number of nonzero for the current column of R and Q
	ScalarVector tval(m);		  // The dense vector used to compute the current column
	RealScalar pivotThreshold = m_threshold;

	m_R.setZero();
	m_Q.setZero();
	m_pmat = mat;
	if (!m_isEtreeOk) {
		m_outputPerm_c = m_perm_c.inverse();
		internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
		m_isEtreeOk = true;
	}

	m_pmat.uncompress(); // To have the innerNonZeroPtr allocated

	// Apply the fill-in reducing permutation lazily:
	{
		// If the input is row major, copy the original column indices,
		// otherwise directly use the input matrix
		//
		IndexVector originalOuterIndicesCpy;
		const StorageIndex* originalOuterIndices = mat.outerIndexPtr();
		if (MatrixType::IsRowMajor) {
			originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(), n + 1);
			originalOuterIndices = originalOuterIndicesCpy.data();
		}

		for (int i = 0; i < n; i++) {
			Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
			m_pmat.outerIndexPtr()[p] = originalOuterIndices[i];
			m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i + 1] - originalOuterIndices[i];
		}
	}

	/* Compute the default threshold as in MatLab, see:
	 * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
	 * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3
	 */
	if (m_useDefaultThreshold) {
		RealScalar max2Norm = 0.0;
		for (int j = 0; j < n; j++)
			max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm());
		if (max2Norm == RealScalar(0))
			max2Norm = RealScalar(1);
		pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();
	}

	// Initialize the numerical permutation
	m_pivotperm.setIdentity(n);

	StorageIndex nonzeroCol = 0; // Record the number of valid pivots
	m_Q.startVec(0);

	// Left looking rank-revealing QR factorization: compute a column of R and Q at a time
	for (StorageIndex col = 0; col < n; ++col) {
		mark.setConstant(-1);
		m_R.startVec(col);
		mark(nonzeroCol) = col;
		Qidx(0) = nonzeroCol;
		nzcolR = 0;
		nzcolQ = 1;
		bool found_diag = nonzeroCol >= m;
		tval.setZero();

		// Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e.,
		// all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at
		// node k. Note: if the diagonal entry does not exist, then its contribution must be explicitly added, thus the
		// trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been
		// found.
		for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp) {
			StorageIndex curIdx = nonzeroCol;
			if (itp)
				curIdx = StorageIndex(itp.row());
			if (curIdx == nonzeroCol)
				found_diag = true;

			// Get the nonzeros indexes of the current column of R
			StorageIndex st = m_firstRowElt(curIdx); // The traversal of the etree starts here
			if (st < 0) {
				m_lastError = "Empty row found during numerical factorization";
				m_info = InvalidInput;
				return;
			}

			// Traverse the etree
			Index bi = nzcolR;
			for (; mark(st) != col; st = m_etree(st)) {
				Ridx(nzcolR) = st; // Add this row to the list,
				mark(st) = col;	   // and mark this row as visited
				nzcolR++;
			}

			// Reverse the list to get the topological ordering
			Index nt = nzcolR - bi;
			for (Index i = 0; i < nt / 2; i++)
				std::swap(Ridx(bi + i), Ridx(nzcolR - i - 1));

			// Copy the current (curIdx,pcol) value of the input matrix
			if (itp)
				tval(curIdx) = itp.value();
			else
				tval(curIdx) = Scalar(0);

			// Compute the pattern of Q(:,k)
			if (curIdx > nonzeroCol && mark(curIdx) != col) {
				Qidx(nzcolQ) = curIdx; // Add this row to the pattern of Q,
				mark(curIdx) = col;	   // and mark it as visited
				nzcolQ++;
			}
		}

		// Browse all the indexes of R(:,col) in reverse order
		for (Index i = nzcolR - 1; i >= 0; i--) {
			Index curIdx = Ridx(i);

			// Apply the curIdx-th householder vector to the current column (temporarily stored into tval)
			Scalar tdot(0);

			// First compute q' * tval
			tdot = m_Q.col(curIdx).dot(tval);

			tdot *= m_hcoeffs(curIdx);

			// Then update tval = tval - q * tau
			// FIXME: tval -= tdot * m_Q.col(curIdx) should amount to the same (need to check/add support for efficient
			// "dense ?= sparse")
			for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
				tval(itq.row()) -= itq.value() * tdot;

			// Detect fill-in for the current column of Q
			if (m_etree(Ridx(i)) == nonzeroCol) {
				for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) {
					StorageIndex iQ = StorageIndex(itq.row());
					if (mark(iQ) != col) {
						Qidx(nzcolQ++) = iQ; // Add this row to the pattern of Q,
						mark(iQ) = col;		 // and mark it as visited
					}
				}
			}
		} // End update current column

		Scalar tau = RealScalar(0);
		RealScalar beta = 0;

		if (nonzeroCol < diagSize) {
			// Compute the Householder reflection that eliminate the current column
			// FIXME this step should call the Householder module.
			Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);

			// First, the squared norm of Q((col+1):m, col)
			RealScalar sqrNorm = 0.;
			for (Index itq = 1; itq < nzcolQ; ++itq)
				sqrNorm += numext::abs2(tval(Qidx(itq)));
			if (sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0)) {
				beta = numext::real(c0);
				tval(Qidx(0)) = 1;
			} else {
				using std::sqrt;
				beta = sqrt(numext::abs2(c0) + sqrNorm);
				if (numext::real(c0) >= RealScalar(0))
					beta = -beta;
				tval(Qidx(0)) = 1;
				for (Index itq = 1; itq < nzcolQ; ++itq)
					tval(Qidx(itq)) /= (c0 - beta);
				tau = numext::conj((beta - c0) / beta);
			}
		}

		// Insert values in R
		for (Index i = nzcolR - 1; i >= 0; i--) {
			Index curIdx = Ridx(i);
			if (curIdx < nonzeroCol) {
				m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
				tval(curIdx) = Scalar(0.);
			}
		}

		if (nonzeroCol < diagSize && abs(beta) >= pivotThreshold) {
			m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
			// The householder coefficient
			m_hcoeffs(nonzeroCol) = tau;
			// Record the householder reflections
			for (Index itq = 0; itq < nzcolQ; ++itq) {
				Index iQ = Qidx(itq);
				m_Q.insertBackByOuterInnerUnordered(nonzeroCol, iQ) = tval(iQ);
				tval(iQ) = Scalar(0.);
			}
			nonzeroCol++;
			if (nonzeroCol < diagSize)
				m_Q.startVec(nonzeroCol);
		} else {
			// Zero pivot found: move implicitly this column to the end
			for (Index j = nonzeroCol; j < n - 1; j++)
				std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j + 1]);

			// Recompute the column elimination tree
			internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
			m_isEtreeOk = false;
		}
	}

	m_hcoeffs.tail(diagSize - nonzeroCol).setZero();

	// Finalize the column pointers of the sparse matrices R and Q
	m_Q.finalize();
	m_Q.makeCompressed();
	m_R.finalize();
	m_R.makeCompressed();
	m_isQSorted = false;

	m_nonzeropivots = nonzeroCol;

	if (nonzeroCol < n) {
		// Permute the triangular factor to put the 'dead' columns to the end
		QRMatrixType tempR(m_R);
		m_R = tempR * m_pivotperm;

		// Update the column permutation
		m_outputPerm_c = m_outputPerm_c * m_pivotperm;
	}

	m_isInitialized = true;
	m_factorizationIsok = true;
	m_info = Success;
}

template<typename SparseQRType, typename Derived>
struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived>>
{
	typedef typename SparseQRType::QRMatrixType MatrixType;
	typedef typename SparseQRType::Scalar Scalar;
	// Get the references
	SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose)
		: m_qr(qr)
		, m_other(other)
		, m_transpose(transpose)
	{
	}
	inline Index rows() const { return m_qr.matrixQ().rows(); }
	inline Index cols() const { return m_other.cols(); }

	// Assign to a vector
	template<typename DesType>
	void evalTo(DesType& res) const
	{
		Index m = m_qr.rows();
		Index n = m_qr.cols();
		Index diagSize = (std::min)(m, n);
		res = m_other;
		if (m_transpose) {
			eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
			// Compute res = Q' * other column by column
			for (Index j = 0; j < res.cols(); j++) {
				for (Index k = 0; k < diagSize; k++) {
					Scalar tau = Scalar(0);
					tau = m_qr.m_Q.col(k).dot(res.col(j));
					if (tau == Scalar(0))
						continue;
					tau = tau * m_qr.m_hcoeffs(k);
					res.col(j) -= tau * m_qr.m_Q.col(k);
				}
			}
		} else {
			eigen_assert(m_qr.matrixQ().cols() == m_other.rows() && "Non conforming object sizes");

			res.conservativeResize(rows(), cols());

			// Compute res = Q * other column by column
			for (Index j = 0; j < res.cols(); j++) {
				Index start_k = internal::is_identity<Derived>::value ? numext::mini(j, diagSize - 1) : diagSize - 1;
				for (Index k = start_k; k >= 0; k--) {
					Scalar tau = Scalar(0);
					tau = m_qr.m_Q.col(k).dot(res.col(j));
					if (tau == Scalar(0))
						continue;
					tau = tau * numext::conj(m_qr.m_hcoeffs(k));
					res.col(j) -= tau * m_qr.m_Q.col(k);
				}
			}
		}
	}

	const SparseQRType& m_qr;
	const Derived& m_other;
	bool m_transpose; // TODO this actually means adjoint
};

template<typename SparseQRType>
struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<SparseQRType>>
{
	typedef typename SparseQRType::Scalar Scalar;
	typedef Matrix<Scalar, Dynamic, Dynamic> DenseMatrix;
	enum
	{
		RowsAtCompileTime = Dynamic,
		ColsAtCompileTime = Dynamic
	};
	explicit SparseQRMatrixQReturnType(const SparseQRType& qr)
		: m_qr(qr)
	{
	}
	template<typename Derived>
	SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other)
	{
		return SparseQR_QProduct<SparseQRType, Derived>(m_qr, other.derived(), false);
	}
	// To use for operations with the adjoint of Q
	SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const
	{
		return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
	}
	inline Index rows() const { return m_qr.rows(); }
	inline Index cols() const { return m_qr.rows(); }
	// To use for operations with the transpose of Q FIXME this is the same as adjoint at the moment
	SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const
	{
		return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
	}
	const SparseQRType& m_qr;
};

// TODO this actually represents the adjoint of Q
template<typename SparseQRType>
struct SparseQRMatrixQTransposeReturnType
{
	explicit SparseQRMatrixQTransposeReturnType(const SparseQRType& qr)
		: m_qr(qr)
	{
	}
	template<typename Derived>
	SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other)
	{
		return SparseQR_QProduct<SparseQRType, Derived>(m_qr, other.derived(), true);
	}
	const SparseQRType& m_qr;
};

namespace internal {

template<typename SparseQRType>
struct evaluator_traits<SparseQRMatrixQReturnType<SparseQRType>>
{
	typedef typename SparseQRType::MatrixType MatrixType;
	typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
	typedef SparseShape Shape;
};

template<typename DstXprType, typename SparseQRType>
struct Assignment<DstXprType,
				  SparseQRMatrixQReturnType<SparseQRType>,
				  internal::assign_op<typename DstXprType::Scalar, typename DstXprType::Scalar>,
				  Sparse2Sparse>
{
	typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
	typedef typename DstXprType::Scalar Scalar;
	typedef typename DstXprType::StorageIndex StorageIndex;
	static void run(DstXprType& dst, const SrcXprType& src, const internal::assign_op<Scalar, Scalar>& /*func*/)
	{
		typename DstXprType::PlainObject idMat(src.rows(), src.cols());
		idMat.setIdentity();
		// Sort the sparse householder reflectors if needed
		const_cast<SparseQRType*>(&src.m_qr)->_sort_matrix_Q();
		dst = SparseQR_QProduct<SparseQRType, DstXprType>(src.m_qr, idMat, false);
	}
};

template<typename DstXprType, typename SparseQRType>
struct Assignment<DstXprType,
				  SparseQRMatrixQReturnType<SparseQRType>,
				  internal::assign_op<typename DstXprType::Scalar, typename DstXprType::Scalar>,
				  Sparse2Dense>
{
	typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
	typedef typename DstXprType::Scalar Scalar;
	typedef typename DstXprType::StorageIndex StorageIndex;
	static void run(DstXprType& dst, const SrcXprType& src, const internal::assign_op<Scalar, Scalar>& /*func*/)
	{
		dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows());
	}
};

} // end namespace internal

} // end namespace Eigen

#endif
